

% compute moments

MomentsZ  = ComputeMoments( ParamsZ  , GridsZ  , LZ  , LhZ  , zetaZ  , vZ  , GEZ  ) ;
MomentsQO = ComputeMoments( ParamsQO , GridsQO , LQO , LhQO , zetaQO , vQO , GEQO ) ;


% Local gains

Factor = 100 ;

Wg_Z  = Factor * ( MomentsZ.Wx  ./ MomentsCE.Wx - 1 ) ;
Wg_QO = Factor * ( MomentsQO.Wx ./ MomentsCE.Wx - 1 ) ;

Ug_Z  = Factor * ( MomentsZ.WUx  ./ MomentsCE.WUx - 1 ) ;
Ug_QO = Factor * ( MomentsQO.WUx ./ MomentsCE.WUx - 1 ) ;

Eg_Z  = Factor * ( ( MomentsZ.WUx  .* MomentsZ.WEx  ) ./ ( MomentsCE.WUx .* MomentsCE.WEx ) - MomentsZ.WUx  ./ MomentsCE.WUx ) ;
Eg_QO = Factor * ( ( MomentsQO.WUx .* MomentsQO.WEx ) ./ ( MomentsCE.WUx .* MomentsCE.WEx ) - MomentsQO.WUx ./ MomentsCE.WUx ) ;

Hg_Z  = Factor * ( ( MomentsZ.WUx  .* MomentsZ.WEx  .* MomentsZ.WHx )  ./ ( MomentsCE.WUx .* MomentsCE.WEx .* MomentsCE.WHx ) - ( MomentsZ.WUx  .* MomentsZ.WEx )  ./ ( MomentsCE.WUx .* MomentsCE.WEx ) ) ;
Hg_QO = Factor * ( ( MomentsQO.WUx .* MomentsQO.WEx .* MomentsQO.WHx ) ./ ( MomentsCE.WUx .* MomentsCE.WEx .* MomentsCE.WHx ) - ( MomentsQO.WUx .* MomentsQO.WEx ) ./ ( MomentsCE.WUx .* MomentsCE.WEx ) ) ;

% Common graph spec

% Area colors
ba = 1 ;
oa = 1 ;
ga = 0.6 ;
gga = 0.4 ;

% x axis and fontsize

fontsize = 14  ;

imin = min( find(uC < 0.2))  - 1 ;
ii = (imin:size(uC))' ;

PercPopC = cumsum(popC) ;

x = flipud( 1 - PercPopC(imin:end) ) ; 
x = x / x(end) ;
xlabe = 'Unemployment rate (p.p.)' ; 
xl = [ -0.03  , 1.03 ] ;
if xl(1) > xl(2)
    xl = fliplr( xl ) ;
end

Treated = GridsZ.tax(ii)-1 > 0.01 ;
TreatFlip = flipud(Treated) ;
xTreat = x(TreatFlip) ;


% Xticks
xt = (0:0.2:1)' ;
xtl = cell(size(xt)) ;
for i = 1:size(xt,1)
    [~,ix] = min( ( ( 1 - x ) - xt(i) ).^2 ) ;
    xtl{i} = num2str( 100 * uC(ix) , 2 ) ;
end
xtl{1} = '6.0' ;
xtl{end} = '20.0+' ;

% Number of subplots
nx=1; ny=2;

% Figure with human capital
figure(4)
set(gcf,'Position',[100 100 1000 400])
yl = [-0.05 * Factor ,0.2 * Factor] ;

subplot(nx,ny,1)
h1 = area(x,flipud([ Ug_QO(ii) + Eg_QO(ii) , Hg_QO(ii) ])) ;
hold on
h2 = area(x,flipud([ Ug_QO(ii) + Eg_QO(ii) , Hg_QO(ii) ])) ;
hold on
ylim(yl)
xlim(xl)
ylabel('Welfare gains (\%)','interpreter','latex')
xlabel(xlabe,'interpreter','latex')
xticks(xt)
xticklabels(xtl)
h1(1).FaceColor = MyBlue ;
h1(2).FaceColor = MyOrange ;
h1(1).FaceAlpha = ba ;
h1(2).FaceAlpha = oa ;
h2(1).FaceColor = MyBlue ;
h2(2).FaceColor = MyOrange ;
h2(1).FaceAlpha = ba ;
h2(2).FaceAlpha = oa ;
title('Quasi-optimal policy','interpreter','latex','FontSize',fontsize)
AxisFonts(0.75*fontsize,fontsize,0.75*fontsize,fontsize)
grid on

subplot(nx,ny,2)
treatedA = area([xTreat(1),xTreat(end)],[yl(2),yl(2)]) ;
treatedA(1).FaceColor = [150 150 150]/255 ;
treatedA(1).FaceAlpha = gga ;
hold on
h1 = area(x,flipud([ Ug_Z(ii) + Eg_Z(ii) , Hg_Z(ii) ]));
hold on
h2 = area(x,flipud([ Ug_Z(ii) + Eg_Z(ii) , Hg_Z(ii) ]));
hold on
ylim(yl)
xlim(xl)
ylabel('Welfare gains (\%)','interpreter','latex')
xlabel(xlabe,'interpreter','latex')
xticks(xt)
xticklabels(xtl)
h1(1).FaceColor = MyBlue ;
h1(2).FaceColor = MyOrange ;
h1(1).FaceAlpha = ba ;
h1(2).FaceAlpha = oa ;
h2(1).FaceColor = MyBlue ;
h2(2).FaceColor = MyOrange ;
h2(1).FaceAlpha = ba ;
h2(2).FaceAlpha = oa ;
title('EZ program','interpreter','latex','FontSize',fontsize)
grid on
legend('EZ locations','Welfare gains given $k=1$','Human capital','Location','NorthWest','Interpreter','latex','FontSize',0.75*fontsize) ;
legend boxoff
AxisFonts(0.75*fontsize,fontsize,0.75*fontsize,fontsize)

if options.SaveResults
    
    % pdf
    filename = 'results/FigureWelfareGainsLocationHumanCap' ;
    run printPDF.m ;
    
    % tiff
    exportgraphics(gcf,'results/FigureWelfareGainsLocationHumanCap.tiff','Resolution',1000) 

end

hold off



% Welfare gains in actual cities
data = readtable('source/LaborMarketBasics.csv');
data.Properties.VariableNames = {'City','EU','UE','W','wt','u'} ;
data.FpSQ = - log( 1 - data.EU - data.UE ) ;
data.s = data.FpSQ ./ ( data.EU + data.UE ) .* data.EU ;
data.fr = data.FpSQ ./ ( data.EU + data.UE ) .* data.UE ;

xF = nan( size(data.s)) ;
yF = nan( size(data.s)) ;
fracF = nan( size(data.s)) ;

for i=1:size(data.s,1)
    [~,xF(i),yF(i),fracF(i)] = FindZeroInterp( data.u(i) - uC , Grids.x ) ;
end

uF = nan( size(data.s)) ;
u0 = nan( size(data.s)) ;

for i=1:size(data.s,1)
    
    u0(i) = (1-fracF(i)) * uC(xF(i)) + fracF(i) * uC(yF(i)) ;
    uF(i) = (1-fracF(i)) * uQO(xF(i)) + fracF(i) * uQO(yF(i)) ...
          - u0(i) ;
end


popF = data.wt / sum(data.wt) ;
WelfareGainsF = 100 * ( MomentsQO.Wx(xF) ./ MomentsCE.Wx(xF) - 1 ) ;

EC = Params.kappa / ( Params.kappa - 1 ) * zetaC(xF) ./ ( zetaC(xF) - 1 ) ;
EP = Params.kappa / ( Params.kappa - 1 ) * zetaQO(xF) ./ ( zetaQO(xF) - 1 ) ;
TFPF = 100 * ( ( ( Params.b + vQO(xF) ) .* EP ) ./ ( ( Params.b + vC(xF) ) .* EC ) - 1 ) ;

if options.SaveResults

    csvwrite_with_headers(...
        '/Users/adrien/Dropbox/JMP/7_statarev/source/WelfareGainsMap.csv',...
        [data.City,popF,u0,uF,WelfareGainsF, TFPF],...
        {'ze','pop','u0','du','welfare','tfp'}) ;

end



